#!/bin/bash
set -e

# while getopts  ":s:p:" opts
# do
        # case  $opts  in
        # s)
		# barcode=$OPTARG
		# sample_name=$OPTARG
		
		# ;;
		# p)
		# out_prefix=$OPTARG
		
		# ;;
        # esac
# done
# shift $(($OPTIND - 1))



if [ -z "$2" ]; then
	echo `basename $0` "<input.vcf> <out.txt>"
	exit 1
fi


#-----------------------------------------------
. /mnt/ilustre/app/medical/tools/.var #---------
#-----------------------------------------------


dat=`date +'%b_%d_%H_%M_%S_%Y'`

fil=`basename $1`

dir_name=${dat}__$fil

mkdir $dir_name
echo $dir_name maked 
echo tep files will be located in $dir_name ...
cd $dir_name


script_path=$tools_path/script
snpeff_path=${tools_path}/snpEff
snpsift=${snpeff_path}/SnpSift.jar

java_memory=1g
echo java memory: $java_memory


# echo
# echo extractField ../$1 CHROM POS
# java -Xmx$java_memory -jar $snpsift extractField ../$1 CHROM POS > 2.txt
echo
echo extractField 
echo

echo CHROM POS
java -jar $snpsift extractField ../$1 CHROM > chrom.txt

echo POS
java -jar $snpsift extractField ../$1 POS > pos.txt

# 将 ID 域预处理，然后将其分离为单独的列
echo
echo id
java -jar $snpsift extractField ../$1 ID | sed -e 's/^COSM/;COSM/' | sed -e 's/\(^rs[0-9]*$\)/\1;/' > 1.id.sep.txt

echo RS ID
cut -d ";" -f 1 1.id.sep.txt > rs_id.txt
echo COSMIC ID
cut -d ";" -f 2 1.id.sep.txt > cosmic_id.txt

echo REF
java -Xmx$java_memory -jar $snpsift extractField ../$1 REF > 1.ref.txt
echo ALT
java -Xmx$java_memory -jar $snpsift extractField ../$1 ALT > 1.alt.txt
cp 1.ref.txt ref.txt
cp 1.alt.txt alt.txt


# paste 2.txt 1.rs.id.txt 1.cosmic.id.txt > 4.txt

# DP REF ALT 需要计算或需要再格式化，所以都要单独提取
# java -Xmx$java_memory -jar $snpsift extractField ../$1 "GEN[0].DP" "DP" | sed -e 's/\t/\//' > 1.dp2.txt
echo DP
java -Xmx$java_memory -jar $snpsift extractField ../$1 "GEN[0].DP" > 1.dp.txt
cp 1.dp.txt dp.txt




# 取得 ref alt 同位基因的深度值 such as 23,22,12
# then merge like this ref=#dp;alt=#dp
echo
echo extractField ../$1 GEN[0].AD[*]
java -Xmx$java_memory -jar $snpsift extractField ../$1 "GEN[0].AD[*]" > 1.ref.alt.dp.txt
paste 1.ref.alt.dp.txt 1.dp.txt > 1.all_dp.txt
java -Xmx$java_memory -jar $snpsift extractField ../$1 "GEN[0].AD[0]" > 1.ref.dp.txt
# -s 保存文件时的域分隔符 -d 从第二位开始取，即为 所有的alt
java -Xmx$java_memory -jar $snpsift extractField ../$1 -s "," "GEN[0].AD[*]" | cut -d "," -f 2- > 1.alt.dp.txt

# format like this: ref=#dp alt1,alt2...,atln=#dp1,#dp2...,#dpn
# then merge like this: ref;alt
paste -d "=" 1.ref.txt  1.ref.dp.txt > 1.ref.dp.merge.txt
paste -d "=" 1.alt.txt  1.alt.dp.txt > 1.alt.dp.merge.txt
paste -d ";" 1.ref.dp.merge.txt 1.alt.dp.merge.txt > 1.all.dp.txt

echo ALLELE COVERAGE
cp 1.all.dp.txt allele_cvg.txt




# 删除第一行的注释用于计算
sed -i "1d" 1.all_dp.txt
echo
echo format depth
perl $script_path/format_0318.pl 1.all_dp.txt > 1.all.frq.txt

# 在第一行插入注释
sed -i "1i depth" 1.all.frq.txt

# 取得 ref 的频率并 format like this: ref=#ref_frq
cut -f 1 1.all.frq.txt > 1.ref.frq.txt
paste -d "=" 1.ref.txt 1.ref.frq.txt > 1.ref_frq.txt

# 取得 alt 的频率并 format like this: alt1,alt2=#alt1_frq,alt2_frq
cut -f 2 1.all.frq.txt > 1.alt1.frq.txt

echo ALT FREQUENCY
cp 1.alt1.frq.txt alt_freq.txt

cut -f 2- 1.all.frq.txt | sed -e 's/\t/,/' > 1.alt.frq.txt
paste -d "=" 1.alt.txt 1.alt.frq.txt > 1.alt_frq.txt
paste -d ";" 1.ref_frq.txt 1.alt_frq.txt > 1.all_frq.txt

# 计算 p-value ??
echo
echo extractField ../$1 QUAL
java -Xmx$java_memory -jar $snpsift extractField ../$1 QUAL > 1.qual.txt
sed -i "1d" 1.qual.txt
echo
echo format p-value
perl $script_path/format_0318.pl 1.qual.txt > 1.p.txt
sed -i "1i QUAL" 1.qual.txt
sed -i "1i P_value" 1.p.txt

echo P-VALUE
cp 1.p.txt p_value.txt

# echo
# echo extractField ../$1 ANN[0].GENE VARTYPE ANN[0].TRID ANN[0].RANK "ANN[0].HGVS_C" "ANN[0].HGVS_P"
# java -Xmx$java_memory -jar $snpsift extractField ../$1 "ANN[0].GENE" VARTYPE "ANN[0].TRID" "ANN[0].RANK" "ANN[0].HGVS_C" "ANN[0].HGVS_P" > i.txt

echo
echo VARTYPE
java -Xmx$java_memory -jar $snpsift extractField ../$1 VARTYPE > vartype.txt

# cut -f1-2 i.txt >i1.txt
# cut -f3 i.txt | enst_refseq.pl -|awk '{print $1 "/" $2}' > i2_trid.txt
# cut -f4- i.txt > i3.txt
# paste i1.txt i2_trid.txt i3.txt > i.new.txt

format_vep.sh ../$1 1.csq.txt
cat 1.csq.txt|cut -d'|' -f2 > consequence.txt
cat 1.csq.txt|cut -d'|' -f3 > impact.txt
cat 1.csq.txt|cut -d'|' -f4 > gene.txt
cat 1.csq.txt|cut -d'|' -f7 > refseq.txt
cat 1.csq.txt|cut -d'|' -f8 > biotype.txt
cat 1.csq.txt|cut -d'|' -f9,10 |sed -e 's/|//' > exon.intron.txt
cat 1.csq.txt|cut -d'|' -f11 > hgvsc.txt
cat 1.csq.txt|cut -d'|' -f12 > hgvsp.txt



# paste 4.txt 1.ref.txt 1.alt.txt i.new.txt 1.dp2.txt 1.all.dp.txt 1.all_frq.txt 1.alt1.frq.txt 1.p.txt > ii.txt


# echo
# echo extractField ../$1 ANN[0].EFFECT ANN[0].IMPACT FILTER dbNSFP_SIFT_score dbNSFP_Interpro_domain dbNSFP_Uniprot_acc dbNSFP_Polyphen2_HVAR_score
# java -Xmx$java_memory -jar $snpsift extractField ../$1 "ANN[0].EFFECT" "ANN[0].IMPACT" FILTER dbNSFP_SIFT_score dbNSFP_Interpro_domain dbNSFP_Uniprot_acc dbNSFP_Polyphen2_HVAR_score > iii.txt

echo FILTER
java -Xmx$java_memory -jar $snpsift extractField ../$1 FILTER > filter.txt

echo Index > index.txt
seq 1 $[`wc -l chrom.txt |cut -d' ' -f1`-1] >> index.txt

sed -i '1c\Chrom' chrom.txt
sed -i '1c\GRCh37 Pos' pos.txt
sed -i '1c\dbSNP ID' rs_id.txt
sed -i '1c\COSMIC ID' cosmic_id.txt
sed -i '1c\Ref' ref.txt
sed -i '1c\Var' alt.txt
sed -i '1c\VarType' vartype.txt
sed -i '1c\Coverage' dp.txt
sed -i '1i\Allele Coverage' allele_cvg.txt
# sed -i '1i\Allele Frequency' alt_freq.txt
# echo Allele Coverage |cat - allele_cvg.txt > allele_cvg1.txt
echo Allele Frequency |cat - alt_freq.txt > alt_freq1.txt
# sed -i '1c\P-VALUE' p_value.txt
echo P-VALUE |cat - p_value.txt > P_value1.txt
sed -i '1c\Filter' filter.txt
sed -i '1c\Gene' gene.txt
sed -i '1c\Transcript ID' refseq.txt
sed -i '1c\Exon/Intron Rank' exon.intron.txt
sed -i '1c\Consequence' consequence.txt
sed -i '1c\HGVS_Codon Change' hgvsc.txt
sed -i '1c\HGVS_AA Change' hgvsp.txt
sed -i '1c\Impact' impact.txt

paste chrom.txt pos.txt rs_id.txt cosmic_id.txt ref.txt alt.txt vartype.txt dp.txt allele_cvg.txt alt_freq1.txt p_value1.txt filter.txt gene.txt refseq.txt exon.intron.txt consequence.txt hgvsc.txt hgvsp.txt impact.txt> ../$2

# paste ii.txt iii.txt > iv.txt

# sed -i "1d" iv.txt
# cat $data_path/title1.txt iv.txt > ../$2



# for sift and polyphen2 format

# considering grchb37 and hg19 are the same in some extent...
# cat 2.txt |awk '{sub(/^/, "chr");sub(/\t/, ":"); print}' > 1.chr_pos.txt
# paste -d "/" 1.ref.txt 1.alt.txt > 1.ref_alt.txt
# paste -d " " 1.chr_pos.txt 1.ref_alt.txt | awk 'NR>1 {print}' > ../1.pph.txt










